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ABSTRACT. In this paper, we discuss a new iterative method for computing sin p . This function 
was introduced by Lindqvist in connection with the unidimensional nonlinear Dirichlet eigenvalue 
problem for the p-Laplacian. The iterative technique was inspired by the inverse power method in 
finite dimensional linear algebra and is competitive with other methods available in the literature. 
Keywords: p-Laplacian, eigenvalues, eigenfunctions, sin p , inverse power method. 

1 Introduction 

In this paper we present a new method to compute the function sin p , inspired by recent 
work done by the authors in [BEM] . where an iterative algorithm based on the inverse 
power method of linear algebra was introduced for the computation of the first eigenvalue 
and first eigenfunction of the Dirichlet problem for the p-Laplacian in arbitrary domains in 
M. N . The functions sin p , 1 < p < oo, can be thought of as generalizations of the familiar 
trigonometric functions. They arise in the unidimensional Dirichlet eigenvalue problem for 
the p-Laplacian and were introduced in this capacity in |Lindqvist] , where a power series 
formula for computing them was also formally given. 

In |BRlj snip functions were utilized to introduce a generalization of the Priifer transfor- 
mation and thus represent, in two phase-plane coordinates, Sturm-Liouville-type problems 
involving the iV-dimensional radially symmetric p-Laplacian L p u := x l ~ N (x 1 ^' 1 \u'\ p ~ 2 u') , 
0^a<s<6<oo. This approach was numerically implemented in [BR2j for an eigenvalue 
problem involving L p with separated homogeneous boundary conditions. In that paper an 
interpolation table for sin p was obtained by numerically solving an ODE. Also in that paper 
the authors raised the question of finding a fast and accurate algorithm for computing sin p . 

Our method depends on the convergence of a sequence of functions whose definition, as 
in [BEM], is motivated by an extension of the inverse power method of linear algebra for 
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obtaining the first eigenvalue and first eigenfunction of finite dimensional linear operators. 
These functions are recursively defined and can be given in integral form, so that they can 
be obtained by numerical integration. 

More specifically, recall that it suffices to obtain sin p in the interval I p = [0,7r p /2], since 
it is extended to the interval [7r p /2,7r p ] symmetrically with respect to tt p /2 and afterward to 
the whole real line R as an odd, 27r p -periodic function (the definition of sin p as well as the 
precise value of ti p are recalled in Section 2). We define the following sequence of (positive) 
functions {<p n } C C l (I p ). Set O = 1 and 

(Vn+lkn+lT 2 ) = -<Pn\(t>n\ P ' 2 if X G I p , 

n+1 (0) = < +1 (n p /2) = 0. 

We prove that the scaled sequence { tfp — 10 n / ll^nlloo} converges uniformly to sin p in I p . 
The functions <p n can be written in integral form as 

i 

<t>n+l(x)=J I J 4>n{s) P ~ dsj d0, X & I p , 

and, therefore, are readily computed using standard efficient numerical methods for definite 
integrals. 

This paper is organized as follows. In Section 2, we recall the definition and some basic 
properties of sin p which will be used in the sequel. In Section 3, we show how to recursively 
construct a sequence of functions which converge uniformly to sin p . Finally, in Section 4 we 
compare the performance of our method with those of |Lindqvist] and [BR2J. 



2 The function sin p 

For the sake of completeness we recall in this section the definition and some properties 
of the function sin p . The unidimensional Dirichlet eigenvalue problem for the p-Laplacian, 
p > 1, is 

4>p (u 1 )' = —Xipp (u) if a < x < b, , * 

u(a) = u(b) = 0, {} 

where ip p (t) = t \t\ p ~ 2 . 

It is easy to verify that if Ai is the first eigenvalue of 

4>p (V)' = -Xip p (v) if a < x < m := (o\ 
v (a) = v' (m) = 0, 



and Vi is the corresponding positive eigenfunction, then Ai is also the first eigenvalue for (JTJ 
with 

v\ (x) if a ^ x ^ m, 

V\ (a + b — x) if m ^ x ^ b, 



Ui (x) 



being the corresponding positive eigenfunction. Moreover, this function is stricly increasing 
on [a, m), strictly decreasing on (m, b] and has only one maximum point which is reached at 
x = m. Thus, Halloo = U\ (m). 

An expression for Ai is well known and can be obtained by integration (see [Otani] ) as 
follows. First multiply (JTJ by u[ and integrate the resulting equation by parts on [a, x] to 
obtain 

rx rx 

ijj p (u[) u'il® — 4> P (u[) u'[dx = —\\ I ip p (ui) u[dx. 



We have 



i[)p(ui) u[dx 



i la 

ui(a) 
u'^x) 



ip p (s) ds 



u 'A a ) 



ip p (s) ds 



Substituting 



whence 



This means that 



ip p (w-l) u[dx 
([5]) and in ([3]) we obtain 

[K - k wn = -a, 



p) P 



— \u[ (a) 


r 




\u (x] 


r 


\u(a)\ p 


P 




P 


\u' (x) 


r 


\vf (a)\ p 


p 




P 


\ui (x) 


r 


\ui (a)\ p 


. p 




P 



(3) 

(4) 
(5) 

(6) 



0. 



?—L\^ + hL\ Ul \p = c, 
p p 

where C is a constant and p' = pj (p — 1) is the conjugate of p. The value of C can be 
found computing the value of this expression at the maximum point m; choosing u\ such 
that U\ (m) = 1 we find 



C 



P~ 1 , / 



V 



u\ (m)\ p + — \u\ (m)\ p — ' 



V 



Therefore, 



(p-1) K^r + Ailm (x)r = A! 

for all x G [a, 6]. 

On the interval [a, m] we have u' ^ 0, hence we can write 



(7) 



u[ (x) 



Ai 
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for all x G [a,m]. Integrating this equation on (a,m) leads to 



b- a / Ai T l(m) ds f 1 ds 



which gives the expression 
where we set 

i f 1 ds 

Making the change of variable s = yt in the last integral and using the classical Beta 
function B we obtain 

n/p 



/ - / tv \l-t) vdt = -B 1--,- 



sin(7r/p) 



(Here one use the properties B(x,y)B(x + y, 1 — y) = x/x sin (Try) and B(l, z) = 1/z with 
x — 1 — 1/p and y = z = 1/p). 
Therefore, 

2^T(vr/p) 
sm(7r/pj 

and 

2^^T(vr/p) 



Ai 



(6 — a) sin (n/p) , 

When a = and b = n p we denote the function ^/p — \ Ul by sin p . Thus, sin p (0) = = 
sin^ (wp/2), Ai = 1 and from (J7J): 

I • zip i sm p 1 

smJ H ^— = 1. 

i p — 1 

It is clear from this equation that sin^, (0) = 1. 

We remark that u = sin p is also the unique solution of the initial value problem 

| u /|p + _H1_ = 1j u (o) = 0j 
p — l 

which can be used to define this function. 

Alternatively, we can define sin p on the interval [0,7r p /2] as an inverse function. In fact, 
multiplying (JS} by tfp — 1 and using (jUJ) with a = and 6 = n p we obtain 



i>sm p (x) ^ g 

/ — i = = x, for x G [0, 7T p /2] , 



I _ SP 

p-1 



that is, snip = ( 1 where 

C(*) == 



ds 



for z G 



o, </p~^ 



With this definition, we extend sin p to the interval [7r p /2,7r p ] symmetrically with respect 
to -Kp/2 and afterward to the whole real line R as an odd, 27r p -periodic function. We list the 
basic properties of sin p : 



1. sin p (0) = = sin p (tt p ), sin p (n p /2) = Hsin^^ = tfp - 1. 

2. sin p (x) is strictly increasing in [0,n p /2] and strictly decreasing in [7r p /2,7r p ] . 



3. |sin;(x)| = {/l 



sin, 



p-1 



3 A sequence uniformly convergent to sin p 

Let I p = [0, vr p /2] and define the following sequence of functions {4> n } C C 1 (J p ). Set (fio = 1 
and 

(VV (0n+l))' = -^p(^n) ifxe/p, 

n+1 (0) = (tt p /2) = 0. 

In this section, we prove that the scaled sequence { tfp — T<j> n / 1 1 0n 1 1 oo J* converges uniformly 
to sin p in I p . Before proceeding, we recall some basic properties of the ip p functions: 

Proposition 3.1. (Basic properties of ip p ) The following holds: 

1. ip p is continuous, strictly increasing and odd, for each p > 1. 

2. V P (ab) = ip p (a) tp p (b) . 



V'p(q) 

^ p (6) 

4 - Op)" 1 = Vy- 

\t\ p 

5 - Jo ^p( s )^ s = — ■ 

By a straightforward calculation we can find the following recursive integral expression 
for the 0„-functions: 



/X I PTTp/2 \ 

Vy I y V^p (0n (s)) rfs j tZ0. 



(12) 



It is clear from fl!2p that each <p n is positive, increasing on I p and reaches its maximum 
value at x = tt p /2. One can obtain an explicit expression for 1; the second function in the 
sequence: 

<i> 1 (x) = j (J p/ ip P (i)ds ) de 



Note that 



7T p 

r \ 2 



0\d0 



■Kp/2—X 

2 



ip P ' (y) dy 



X 



1101 



7T, 



1 / 7T p \ P p — 1 f 71 /p 



2 J p V 2 J p \sin(7r/p) 

The next n -functions however, are very difficult to obtain explicitly by solving the integrals 
analytically. On the other hand, the integrals can easily be solved numerically. 



Proposition 3.2. n+ i ^ 



'1 no Vn 



on I n . 



Proof. For n = 1 the result is trivially true since 0o = 1. Assuming by induction that 
0n ^ H0i Hoc 0n-i, we have 

rx / pTTp/2 

n+ l (x) = J^P'lj ^ n ( S )) ds I d0 

TTp/2 



Vy I Vv(ll0ilL) / r l ,lo„- i («}) f/.s j r/0 

Il0l|loo0" ( X ) • 



The following technical lemma, which will be used in the sequel, can be proved via the 
Cauchy mean value theorem (see |AV V] ) and works as a L'Hopital's rule in order to get 
monotonicity for a certain quotient function. 
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Lemma 3.3. Let f,g : [a, b] — > R be continuous on [a,b] and differentiable in (a,b). 

, , f(x) - f(a) , f(x) - f(b) , . , . . . r , 

both —— — ana ——- 7TT are (strictly} increasing [decreasing]. 

g(x) - g(a) g{x) - g(b) 

Theorem 3.4. For each n ^ 1 the function -j 



is strictly decreasing on I p and 



(i) 1 < in f - ^ ^p/ 2 ) 



1101 



4> 



'n+l 



I P <p n+ l <p n+ i(n p /2) \\<p n +l 
Io P/2 (<f>n-l (s))ds 



Io P/2 1pp(<f>n(s))ds 



for n > 1. 



(iii) 


4>n 




0n-l 




• < 


01 










4>n+l 


oo 


0n 


oo 




02 



< oo. 



Proof. Since <pi is strictly increasing, it follows that l/(f>i is strictly decreasing. Assume by 
induction that 4>n-i/4>n is strictly decreasing. Since 

(j) n (x) - (f) n (0) (j) n (x) 



n+ l - (f) n+1 (0) <f) n+1 (x) ' 

in order to show that n /0 n+ i is strictly decreasing, it suffices in light of the lemma to verify 
that </4/</4 + i is strictly decreasing on I p . But, 



0n+l (X) 



7Tp/2 



il; p (<f> n -i(s))ds 



( f* p /2 ^ 
Ipp (0„_i (s)) 



7Tp/2 



V 



7Tp/2 



/ 



Since is strictly increasing and the functions JJ P ^ 2 ip p (0 n -i ^ s an d JJ P ^ 2 Vv (0n (s)) 
are null at x — vr p /2, we can apply the lemma again to verify that the quotient of these 
integral functions is a strictly decreasing function. We have 



7Tp/2 



Vv(0n-1 {s))ds 



7Tp/2 



i) p {<Pn{s))ds 



jM0n-l (g)) 
(0n («)) 



which is strictly decreasing by the induction hypothesis. 
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The inequality in (i) follows from Proposition 3.2. Before verifying (ii) we remark that 
II II oo = 00 smce 0i (0) = 0. In order to prove (ii) we first observe that the monotonicity 
of n /0n+i implies that 



(fir. 



0n+l 



(fin (x) 



00 0n+l ( X ) 



L'Hopital's rule then yields 



lim 



^ n ^ - lim ^ n ^ -"'<-t ^ P ^ n_1 ^ ' 



*To+ <fi n+1 (x) x "6+ (x) rP ' ^ J o ^ /2 ^ (0 n ( s )) rf s 
The proof of (iii) is a consequence of the following estimates, valid for n ^ 2: 



< 00. 



0n 







n+1 



(0n-l (s)) tfe 





/ lfi p ((fi n (s)) ds 
./O 



lfip ((fin s )) 



^ {S) ) dS 



\ 

I f-K p /2 



PTVp/2 

/ Ipp {(fin (s)) ^ 
./0 



S$ ^P 



^p(0n(s))^p 



>n-\ 



ds 



\ 



7Tp/2 



>n-\ 



(fin 
0n-l 



^p 



^P (0n (s))tZs 



/ /-Tp/2 \ 
ifip ((fin (s)) ds 





7Tp/2 



(0 n (s)) C?S 
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Theorem 2.4. Let u n : = 



Il0n|loo 

decreasing for each x E I p and 



G C 1 (Ip) , for n ^ 1. T/ien £/ie sequence {u n (x)} n>1 is 



\fp — lu n — > sin p uniformly in I p . 
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Proof. In I p we have 



u, 



U n+ i n+ l \l|0n+l 



> inf 



-i 



'p (Pn+l 



>n+l\ 



-1 



J n+1\ 



-1 



that is, {u n (x)} n>1 is decreasing for each x G I p , and the whole sequence is bounded below 
by u x . Thus, there exists 

u := \imu„. 



We have ||w n |loo = 1 f° r each n. Moreover, since 



110, 



H0n+l|L 

we also have, for every x G I p , 
\<(x)\ 



inf 



4>n 




01 






0n+l 


oo 


02 



=: C, 



-A 



110, 

ll0nlL 



n p /2 



7Tp/2 



V> P (0n-i(s))da 



110 



n-l 



"TTp/2 

^ C'ty ( J A(l)ds 



It follows from Arzela-Ascoli's theorem that u n — > u G C (I p ), uniformly 
In order to conclude the proof, we need just to show that 



sin, 



u 



From (1121) we can write the following expression: 

r.TTp/2 

U n +1 (x) =ln I i>- 



ip p (u n (s)) ds d9, 



(13) 



where 

^ . = Un\\oo 

Il0n+l|loo' 

In view of the boundedness of {7„}, there exists 7 := lim7 nfc for some subsequence {7„ fc } • 
Thus, letting k — > 00 in 

fTTp/2 

u nk+ i (x) = 7„ fc / Vy I / ( U n k 0)) ds I de, 



we get 

u (x) = 7 



io ^ \l e ^^{^dsUeeC 1 ^), 



which means that u is a positive solution to the following problem 

ip P (u')' = -Jip p (u) if x e I p , 
u (0) = u' (tt p /2) = 0. 

In view of the positivity of u, we can integrate the equation above multiplied by vl and 
proceed as in the derivation of (Q to find 7 = 1. From this we conclude that in fact 
lim7„ = 1 (the whole sequence converges to the eingenvalue 1) and that u = limM„ satisfies 
the same boundary value problem that sin p / ||sin p || does. Since both u and sin p are positive 
and = ||sin D / ||sinJ||l = 1, we must have 

II Moo 11 p I \\ p\\ lloo " 



U 



sm p 
I sirip I 



whence (TT~3T) follows. 



4 Numerical Results 

Next we examine the computational time of each method. Computations were performed on 
a WindowsXP/Pentium 4-2.8GHz platform, using the GCC compiler. Although the method 
of computing sin p by solving an ODE suggested in [BR2j (which we implemented by means 
of a standard Runge-Kutta fourth power method) is by far the fastest, the computational 
times of the other two methods are competitive, the inverse power method being on average 
more than twice as fast as the power series method of ILindgvist] for values of p greater than 
2. Also, the average number of 8 iterations that the inverse power method uses to obtain 
the same (and sometimes better; see Table 2) accuracy of the differential equation method 
of [BR2j is quite remarkable, specially taking into account that the functions (f) n converge to 
rather rapidly. We emphasize that the computational time of the inverse power method is 
not the main subject of this presentation. The method demands the computation of double 
integrals at each iteration for each grid point. We opted for a classical, computationally easy 



10 



to implement and reasonably fast method to compute these integrals, namely, the Simpson 
composite method. However, a greater effort spent in lessening the computational time of 
the numerical integrations certainly would be reflected in a substantial decrease in the time 
spent computing sin p overall. Nevertheless, by considering the accuracy and the comparison 
scale among the three methods (on the range of miliseconds) we may say that the results 
presented in this paper validate the inverse power method as an effective and reasonably fast 
method for numerically obtaining sin p . 

Below we present the average time spent in computing sin p on the whole interval I p 
divided in 101 grid points by each method for six values of p (the average was taken out 
of five computer runs); the stop criterion in each method was an error tolerance of 1CT 8 
between successive iterations and less than 500 terms in the power series. 



p 


1.1 


1.5 


2.0 


2.5 


3.0 


3.5 


Inverse power method 


21.5 


32.1 


1.1 


37.7 


37.8 


31.7 


Differential equation method 


1.9 


1.8 


1.1 


1.5 


1.5 


1.5 


Power series 


92.9 


2.2 


2.0 


79.6 


79.3 


73.3 



Table 1: Average time (in miliseconds) for the computation of sin p on I p for each method. 

Besides the trivial point 0, the only point where the value of sin p is exactly known is 7r p /2, 
with sin p (vr p /2) = tfp — 1. In the next table we present the computed value for sin p (7r p /2) 
obtained using each method: 



p 


1.1 


1.5 


2.0 


2.5 


3.0 


3.5 


</p-l 


0.123285 


0.629961 


1 


1.17608 


1.25992 


1.29926 


Inverse power method 


0.123285 


0.629961 


1 


1.17608 


1.25992 


1.29926 


Differential equation method 


0.123285 


0.629966 


1.00017 


1.17647 


1.26044 


1.29983 


Power series 


5.3 x 10 128 


0.629961 


1 


1.17608 


1.25993 


1.29928 



Table 2: Value of sin p (tt p /2) = tfp — 1 obtained independently using each method. 



Notice that the inverse power method appears to be more accurate when computing sin p at 
values close to 7r p /2. Indeed, in order to obtain a good approximation close to this point, 
it was necessary to allow for a greater number of terms in the power series than would be 
necessary for points far from 7r p /2. 



p 


1.1 


1.5 


2.0 


2.5 


3.0 


3.5 


Inverse power method 


5 


8 


9 


8 


8 


8 


Power Series 


501 


13 


8 


470 


501 


501 



Table 3: Number of iterations. 
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We see that the number of iterations used by the inverse power method is remarkably low. 
Below, we present the graphics of sin p for the same values of p computed using the three 
methods (except for p — 1.1, since the power series appears to diverge in this case). Notice 
that all three methods agree very well with each other, being virtually indistinguishable. 




0.5 1 1.5 2 2.5 0.5 1 1.5 2 2.5 3 3.5 
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